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Motivated by recent experiments on two-component systems, we investigate the ground-state 
phase diagram of a mixture of two bosonic species by means of path-integral quantum Monte Carlo 
by the two-worm algorithm. The mixture is trapped in a square lattice at different filling con¬ 
ditions. Various quantum phases are stabilized depending on the interplay between intra- and 
inter-species interactions and on the filling factors. We show that the ground-state phase diagram 
at half-filling features a demixed superfluid phase and demixed Mott-Insulator phase when the 
inter-species interaction becomes greater than the intra-species repulsion, and a double-superfluid 
phase or a supercounterflow otherwise. We show that demixing, characterized by spatial separa¬ 
tion of the two species, can be detected experimentally through the effects of anisotropy revealed 
by time-of-flight images. We also study how demixing effects depend on the filling factor of the 
two components. Finally, we found that super-counterflow phase is preserved in the presence of 
unbalanced populations. 


I. INTRODUCTION 

Mixtures of two bosonic species trapped in optical lat¬ 
tices feature a variety of unprecedented effects and quan¬ 
tum phases [l]-[25] resulting from the interplay between 
kinetic energy and intra-, inter-species density-density in¬ 
teractions. In the past decade, a considerable theoretical 
work has been devoted to investigating manifold proper¬ 
ties of these systems. Different aspects of the phase di¬ 
agram have been studied by means of generalized mean- 
field schemes [15]-[17], Luttinger-liquid picture [18], or 
perturbation methods [19]. Moreover, the effect of phase 
separation [20, 21], the study of quantum emulsions and 
coherence properties of mixtures [22, 23], and the shift 
of Mott domains due to the presence of a second (super¬ 
fluid) species [24] along with the interpretation of this 
shift in terms of polaron excitations [25] have been ex¬ 
plored. 

In the limit of large interactions, magnetic-like phases 
such as the incompressible double-checkerboard solid and 
the supercounterflow have been predicted theoretically 
[l]-[3] at total Ailing one and repulsive inter-species in¬ 
teraction, while paired superfluidity has been found [4] in 
the case of attractive inter-species interaction at equal in¬ 
teger filling of the two components. These findings stim¬ 
ulated further investigation of magnetic-like phases in¬ 
cluding finite temperature effects [9, 13], different optical 
lattice geometries [14] and dimensionality, and various 
interacting regimes [5]-[12], [19]. Nevertheless, over ten 
years from initial theoretical investigation of these sys¬ 
tems [l]-[3], their rich phase diagram still exhibits several 
unexplored aspects that challenge theoretical and numer¬ 
ical techniques, while its elusive character demands more 
sophisticated experimental techniques for the observation 
of these quantum phases. 


The recent experimental realization of mixtures, either 
combining two different atomic species [26, 27] or using 
the same atomic species in two different internal energy 
states [28, 29, 37, 38] demonstrated how refined experi¬ 
mental techniques allow to control the model parameters, 
hence reinforcing the interest toward these systems. The 
possibility to observe such new phases in real systems is 
strongly affected by i) the presence of the trapping poten¬ 
tial which introduces an undesired source of inhomogene¬ 
ity, ii) the fact that their theoretic prediction is based on 
assuming rather ideal conditions (such as, for example, 
species A and B with the same boson numbers Na = N^, 
or large intra-species interactions), and Hi) the difficulty 
in reaching low enough temperatures where such phases 
are expected. Concerning points i) and ii), the experi¬ 
mental realizability of magnetic phases has been analyzed 
in [30] leading to promising results at least for the dou¬ 
ble checkerboard phase. In Ref. [30] it was shown that 
the double checkerboard state can be found for large but 
finite intra-species interaction and for certain parabolic 
confinements and particle number imbalance. 

In this work, we reconstruct the ground-state phase 
diagram of a mixture of two twin species (two bosonic 
components with equal hopping parameters and equal 
intra-species interactions) at half-filling ^ and deter¬ 
mine under which conditions the supercounterflow and 
demixed phases are stabilized. We further explore demix¬ 
ing away from v = ^. Demixed phases are characterized 
by spatial separation of the two species. So far, demixing 
effects have been mainly studied in the context of con¬ 
tinuous systems [31, 39-43, 45] and Bose-Fermi mixtures 
[46, 47]. Here, we study demixing in a binary mixture of 
bosons described by the two-component Bose-Hubbard 
model. At ^ we find a demixed superfluid (dSF) 
phase or a demixed Mott-Insulator (dMI) phase when 
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the inter-species interaction is greater than the intra¬ 
species repulsion, and a double-superfluid (2SF) phase 
or a supercounterflow (SCF) otherwise. We characterize 
transitions to demixed phases by introducing a suitable 
demixing parameter. Further significant information is 
found by looking at the off-diagonal correlator which al¬ 
lows to design experimental observation of the various 
phases through time-of-flight images. We also study SCF 
away from = | and find that SCF is stabilized for com¬ 
mensurate total filling although the superfluid response 
in the counter-flow channel does depend on the popula¬ 
tion imbalance. This paper is organized as follows. In 
section II we present our model and formalism. In sec¬ 
tion III we discuss the ground state phase diagram at 
v = supercounterflow properties for unbalanced popu¬ 
lations, and demixing effect at non-commensurate filling. 
Finally, in section IV we conclude. 

II. METHOD AND MODEL 

We study a mixture of two bosonic species in a uniform 
two-dimensional (2D) square optical lattice. The system 
is described by the two-component Bose-Hubbard (BH) 
model: 

^ ~ “f “f ^ab ^ ^ '^ai'^bi (1) 

i 

where Uab is the inter-species repulsion, Uai^bi is the num¬ 
ber operator at site i for species A and B, and 

He = ^ ^ blciijlci 1) tc ^ ^ •) (^) 

* (b> 

with c = a,b denoting the bosonic species, and operators 
Ci, c| satisfying the standard commutator [ci,cj] = 1. 
Parameters Uc represents the intra-species repulsion and 
tc is the hopping amplitude for component c. The symbol 
(ij) refers to summation over nearest neighboring sites. 
To further simplify the number of free parameters, we 
work with twin species, that is, we set Ua = Ub = U, 
ta = tb = t- We are interested in exploring the phase 
diagram of model (1) as a function of U/t and Uab/t 
with particular emphasis on the supercounterflow and 
demixed phases. Our results are based on large-scale 
path-integral quantum Monte Carlo simulations by a 
two-worm algorithm [8]. Unless otherwise noted, we per¬ 
form simulations for system sizes L = 8, 16, 24, 36 (we 
use the lattice step A/2 as a unit length) and we work 
at inverse temperature (5 = L/t which ensures that the 
system is in its ground state. 

III. RESULTS 

Ground-state phase diagram at half-filling. The 
ground-state phase diagram of model (1) is shown in 
Fig. 1 in the U/t vs Uab/t plane. The 2SF phase fea- 



FIG. 1. (Colors online) Ground-state phase diagram of 
twin bosonic species trapped in a uniform optical lattice at 
i/ = I. Four phases are stabilized: double superfluid (2SF), 
supercounterflow (SCF), demixed Mott-Insulator (dMI), and 
demixed superfluid (dSF). Symbols mark the phase bound¬ 
aries as calculated from Monte Carlo simulations (see text). 
Whenever not visible, error bars lies within the symbol size. 
The inset shows the scaled total superfluidity ptot as a func¬ 
tion of the inter-species interaction for system sizes L = 
8,16,24,36 and U/t = 20. The intersection between curves 
corresponding to different system sizes gives the transition 
point for the 2SF-SCF transition. 


tures two 17(1) broken symmetries and is characterized 
by order parameters (a) 0 and (b) 0, or, equivalently, 

finite stiffness of the total superfluid flow ptot 0, and 
finite stiffness of the relative superffuid flow psCF y^ 0. 
The total and relative superfluid stiffnesses are given by 
Ptot.sCF = {{Wa±Wbfi)/2j3, where Wa,b are winding 
numbers of worldlines of species A and B [50]. The SCF 
phase restores one 17(1) broken symmetry and is charac¬ 
terized by order parameter {aU) 0 while (a) = 0 and 
(6) = 0, or, equivalently, zero total superfluid stiffness, 
Ptot = 0, and finite relative superfluid stiffness, psCF 0. 
The demixed phases are characterized by spatial separa¬ 
tion of the two components. This phenomenon is ob¬ 
served whenever Uab > U (as found in [1] within the 
isospin picture of bosonic mixtures for the Mott region). 
A heuristic derivation of this condition for the case of 
generic filling factor is given in the Appendix A. In a 
demixed phase the density distribution of the two com¬ 
ponents is anisotropic on the lattice with density max¬ 
ima of one species corresponding to density minima of 
the other. The dSF features two 17(1) broken symme¬ 
tries, but the two species occupy different regions of the 
lattice. Finally, the dMI is an incompressible, insulating 
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phase where the two components are spatially separated. 
Overall, 2SF and dSF are conducting phases, while dMI 
and SCF are insulating phases, although SCF supports 
flow in the so-called particle-hole channel. In agreement 
with [44], we find that demixing is observed as soon as 
Uab/t > U/t for any value of U. It is also worth noting 
that in order to reach the insulating phases, both intra- 
and inter-species interactions have to be large enough. 
Indeed, for Uab/t < 13.5 {U/t < 15.5) the 2SF (dSF) 
phase is stable for arbitrarily large U/t (Uab/t). 

In order to extract the phase diagram shown in Fig. 1 
we measure superfluid stiffness in terms of winding num¬ 
bers statistics and the demixing parameter A (see below) 
which depends on the density distribution. Both observ¬ 
ables are readily available within the path-integral formu¬ 
lation by the worm algorithm. The demixing parameter 
is given by: 

A_ 1 r (’^ai) — 1 ^ 

- M ^[{na^) + {nb.)\ 


where {nd) is the quantum-thermal average of the den¬ 
sity of component c = a, 6 at site i. Parameter (3) is 
basically a lattice average of the square of the imbalance 
of the local species’ density. This parameter is different 


than the one used in Ref. [31] D = 


(jVa)-(jVb) 
{N^) + {Nt) 


but 


gives the same information about demixed phases. 

In Fig. 1, solid circles correspond to the 2SF-SCF tran¬ 
sition. This transition belong to the (2-|-l)-XY univer¬ 
sality class. Transition points are found using standard 
finite-size scaling of ptot as it can be seen in the in¬ 
set of Fig. 1 where we plot the scaled total superfluid¬ 
ity as a function of interaction Uab/t for system sizes 
L = 8,16,24,36 and U/t = 20. The curves correspond¬ 
ing to different sizes intersect at the critical point Uab/t 
= 14.9 ± 0.1. Both 2SF and SCF are stable for Uab < U. 

We detect the phase transition between 2SF and dSF 
(squares in Fig. 1) by studying the behavior of the A 
parameter. In Fig. 2 we show A as a function of Uab/t 
for U/t = 5, U/t = 10, U/t = 15, circles, squares, trian¬ 
gles respectively. A jump of four orders of magnitude is 
clearly visible when Uab/t ^ U/t signaling the onset of 
the demixed phase dSF at the expense of the 2SF. Fur¬ 
ther significant information about the demixed phase is 
achieved by looking at the density distribution of par¬ 
ticles in the lattice. The quantum-statistical average of 
the particle number (uai), (nbi) is displayed in Fig. 3. 
The X and y axis denote the x and y coordinates on the 
lattice. The color code is displayed on the the right bar. 
Panel a) refers to the 2SF case where the density of each 
species is uniformly distributed in the lattice, correspond¬ 
ing to the spatial coexistence of the two species. Panel 
b) refers to the dSF phase. Here, we clearly see that the 
two components occupy spatially-separated regions with 
well defined boundaries of a few lattice steps of thickness 
where the two components coexist. 

The dSF phase becomes unstable towards a dMI phase 
upon increasing the intra-species interaction (triangles in 



FIG. 2. (Colors online) Demixing parameter A as function 
of Uab/t for U/t = 5, 10,15 (circles, squares, triangles respec¬ 
tively) across the 2SF-dSF transition. A jump of four orders 
of magnitude is clearly visible when Uab/t ~ U/t signaling the 
onset of the demixed phase dSF at the expense of the 2SF. 
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FIG. 3. (Colors online) Quantum-statistical average of the 
particle number (uai), {ribi}. The x and y axis denote the x 
and y coordinates on the lattice. The color code is displayed 
on the the right bar. Panel a): 2SF phase at U/t = 10, 
Uab/t = 7 with density uniformly distributed in the lattice. 
Panel b): dSF phase at U/t = 15, Uab/t = 20 with two com¬ 
ponents occupying spatially-separated regions of the lattice. 
Panel c): dMI phase Sut U/t = 20, Uab/t = 22 with two com¬ 
ponents occupying spatially-separated regions of the lattice. 
Compenetration between the two regions is noticeable. 
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Fig. 1). We compute the phase boundary by verifying 
the drop of the SF density of each component for system 
sizes L = 8, 16, 24. While the dSF is characterized by two 
U(l) broken symmetries over spatially separated regions, 
in the dMI these symmetries are restored, that is, the sys¬ 
tem loses its off-diagonal long range (anisotropic) corre¬ 
lations and becomes insulating. Investigating the details 
of this phase transition is challenging. Finite-size scaling 
of the SF stiffness cannot be performed in proximity of 
demixed regions since the statistics of winding numbers 
is affected by the topography and the (non-connected vs. 
connected) topology of demixed regions, both of which 
depend on the initial conditions. The averaged particle 
number within the lattice in the dMI phase is displayed 
in panel c) of Fig. 3. Unlike the dSF phase (panel b)), 
the boundaries between the regions occupied by the two 
species tend to be rough and irregular, and compenetra- 
tion between the regions is more pronounced due to the 
reduced mobility of bosons in the dMI phase. 

Finally the empty circles in Fig. 1 correspond to the 
SCF-dMI transition. Upon entering the dMI the sys¬ 
tem restores the U(l) broken symmetry characterizing 
the SCF. We measured the A parameter and psCF across 
the transition line and, similarly to the 2SF-dSF tran¬ 
sition, we observe an increase of about three orders of 
magnitude in A, while psCF goes to zero. 

In figure 4 we show the computed momentum distri¬ 
butions rick = |</'(k)P Jhij [51] for species 

c = a, b proportional to Time-Of-Flight (TOF) images 
detectable experimentally. Here (()(k) is the Fourier trans¬ 
form of Wannier function <^(r), which we do not compute 
here. TOF profiles along x and y lattice directions within 
the first Brillouin zone are plotted in figure 4 for the 2SF 
(panel a)), dSF (panel b)), SCF(panel c)) and dMI (panel 
d)). Note that the TOF image of SCF corresponds to 
the one of the particle-hole pair. The insets show the 
corresponding quantum-statistical average of the density 
of the two components within the lattice. We observe 
different profiles featuring the different quantum phases 
with an evident detectable anisotropy of distributions rick 
along the x— and y—direction for dSF. This is expected 
due to the anisotropy of the spatial separation and rep¬ 
resents a clear experimental signature of the dSF. On 
the other hand the TOF image of dMI doesn’t reflect 
anisotropy of spatial separation due to the absence of 
off-diagonal long range order. In general, we observe no 
differences in the momentum distributions between the 
two species. 

Non-Commensurate Filling. We now turn to studying 
model (1) away from commensurate total filling where 
neither the dMI nor the SCF exist. In particular, we are 
interested in studying the properties of the dSF phase as 
a function of the total filling n. We find that demixing 
effects are still present at n < 1 although the spatial sep¬ 
aration between the two components A, B is not as pro¬ 
nounced as for n > 1. This can be seen in Fig. 5, where 
we plot the demixing parameter A as a function of n (cir¬ 
cles). A substantial drop in the value of A is observed 



FIG. 4. Momentum distributions Uk in the first Brillouin 
zone for double superfluid 2SF (panel a)), demixed super¬ 
fluid dSF (panel b)), supercounterflow (panel c)) and demixed 
Mott-Insulator (panel d)). The insets show the corresponding 
quantum-statistical average of the density of the two compo¬ 
nents within the lattice. The demixed phases dSF shows an 
evident anisotropy along the x— and i/—direction. This is 
expected due to the anisotropy of the spatial separation. 


for n < 1. This is due to the presence of large regions in 
the lattice where A and B overlap, as shown in the left 
inset of Fig. 5 where we plot the quantum statistical av¬ 
erage of the densities of the two species at n = 0.3. This 
overlap region results from enhanced hopping of particles 
which is responsible for larger fluctuations of {ua^b)- For 
comparison, in the right inset of Fig. 5, we show the 
species densities at n = 1.3. A net spacial separation 
between A and B is observed. Despite this substantial 
drop in the value of A for n < 1, we find that A is still 
a good indicator that demixing has occurred. Indeed, A 
values in the 2SF phase (triangles in Fig. 5) are still or¬ 
ders of magnitude smaller than in the dSF phase. These 
results suggests that demixing effects can be observed in 
the presence of an external harmonic trap where a vari¬ 
ation of n within the trap is present. See Supplemental 
Material for further details. 


Unbalanced Populations. We conclude by studying the 
SCF in the presence of population imbalance. The SCF 
phase can be stabilized at n = 1. Here we are interested 
in showing that SCF still exists with non-zero imbalance 
although its robustness depends on the latter. In Fig. 
6 we plot pscF for different values of ratio N^/Na. We 
observe that SCF remains robust also for large popula¬ 
tion imbalance although the largest superfluid response 
corresponds to the balanced case. 
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FIG. 5. A parameter as a function of total filling n for dSF 
(circles, U/t = 10, Uab/t = 15 and L — 24) and 2SF (trian¬ 
gles, U/t = 10, Uab/t = 7 and L = 24) phases. Insets repre¬ 
sent the average densities of the two components for n = 0.3 
(left), n = 1.3 (right). 
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FIG. 6. Super-counterfiuidity stiffness as a function of Nb/Na 
{U/t = 20, Uab/t = 17 and L = 24). Dashed line is a cubic 
polynomial least-square fit, while dotted lines represent one-u 
deviation contours from the fit. 


IV. CONCLUSION 

Motivated by recent experiments on two-component 
systems, we have investigated properties of a mixture 
trapped in a 2D optical lattice at different filling fac¬ 
tors. We have used path-integral quantum Monte Carlo 
by the two-worm algorithm. Ai v = i, we have observed 
a demixed superfluid phase and demixed Mott-Insulator 
phase when the inter-species interaction becomes greater 
than the intra- species repulsion, and a double superfluid 
phase or a supercounterflow otherwise. Demixing is char¬ 


acterized by spatial separation of the two components 
and manifests itself experimentally with anisotropic time 
of flight images. The latter have been calculated for se¬ 
lected examples, comparing the results for demixed and 
non-demixed phases. Finally, we have shown that the 
SCF phase survives in the presence of population imbal¬ 
ance and trapping potential (see Supplemental Material). 
In the future, we plan to investigate how finite tempera¬ 
ture affects demixing. Our preliminary results show that, 
in the dSF phase, thermal fluctuations destroy demixing 
effects earlier than they destroy superfluidity, leaving the 
system in a uniform superfluid phase. Moreover, since 
preliminary results show that the demixing parameter is 
sensitive to temperature, this suggests its possible appli¬ 
cation as a thermometer for mixtures of ultracold gases. 
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Appendix A: A simple interpretation of demixing 
effect 

The mixing/demixing effect can be interpreted in a 
simple way for the dSE-SE transition at generic filling fac¬ 
tor. In a perturbative framework, where the energy con¬ 
tribution of tunneling processes is assumed to be negligi¬ 
ble, one should compare the energy of the ground-state 
with the two species spatially separated to the energy 
of the ground-state where the two species coexist within 
the lattice. When the two species occupy spatially sepa¬ 
rated regions of the lattice TZa and TZb, we assume that 
Ma {Mb) sites of TZa {TZb) are occupied by bosons A (B), 
with Ta {rb) sites containing n + \ {m + 1) bosons, and 
Ma — Ta {Mb — Tb) sites containing n {m) bosons. Note 
that the total number of sites is given by M = Ma + Mb 
while the total number of particles is Na = MaU + Ca, 
Nb = MbUi + rb. The nonuniform filling in TZa and TZb re¬ 
flects the SE character of the two species. The resulting 
energy reads: 

Eo= ^Man{n - 1) -f ^Mbm{m - 1) 

+ UaTaU + UbUm - fiaNa - fJ-bNb , (Al) 

where the t/ab-dependent term is absent due to the spa¬ 
tial separation of the two species. The mixing effect is 
described by: a boson is lost from each site of TZa {TZb) 
occupied by n -I- 1 {m + 1) bosons while Va {cb) sites ap¬ 
pear in TZb (TZa) with m bosons B and one boson A {n 
bosons A and one boson B). Uab interaction term is now 
activated and the resulting energy is 

K= ^Man{n - 1) -f ^Mbm{m - 1) 

+ Uab{ram + rbn) - HaNa-fJ-bNb. (A2) 
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This mutual exchange of bosons between TZa and TZb rep¬ 
resents the mixing process with the lowest-energy cost in 
the minimum-energy scenario. The condition Eq < Eq 
(justifying the transition from the uniform ground state 
to the demixed state) implies that UabiTa^n r^n) > 


UaTaiT- + UbThTTi which reduces to the well-known condi¬ 
tion Uab > U for n = m and Ua = Ub- This elementary 
argument is valid in the SF regime due to its semiclassical 
character. It cannot be extended to the transition from 
the dMI-SCF phase where quantum correlations and hop¬ 
ping processes play a preminent role. 
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Here we present further details in our analysis on quantum demixing effect in mixtures 
of two bosonic species. In particular, we wish to provide further information for the case 
of non-commensurate hlling, analyzing the system in non-homogeneous conditions in the 
presence of a harmonic trap. In the following we consider balanced populations. In the 
presence of a harmonic potential, the chemical potential of species c = a,b, occurring in the 
free energy F = H — Yhct^cindi transforms according to 

hci = /ic - (1) 

where uh is the strength of the harmonic trap (expressed in unit of the hopping amplitude 
t), and fi the position vector of lattice site i. This leads to a variable local hlling factor 
Uff through the lattice, allowing us to explore in a more complete way the inhuence of the 
hlling variation on the phase diagram of the system. As we showed in the main paper, 
demixing persists away from hlling n = 1, as long as Uab > U. This is even more evident 
in the trapped case. In hgure 1 we show the density map of a demixed superhuid in which 
a harmonic conhnement forces the mixture to occupy only the central regions of the lattice 
with radially decreasing hlling. In the right panel of hgure 1 we display the demixed density 
prohles of the two species for a radial section perpendicular to the interface between the two 
superhuids. It is evident that demixing persists at any value of the radial coordinate and 
for diherent values of hlling. It is worth noticing that, as predicted in [1] for a mixture of 
two Bose-Einstein condensates in a harmonic trap, the overlap energy is minimized when 
the boundary is a straight line. 

Furthermore, for higher values of the inter- and intra-species potentials, within the SCF 
region of the phase diagram (Fig. 1 in the main text), we observe the presence of the well- 
known spatial shell structure [2-4]. For increasing total number of particles in the system, 
we observe a spatial structure in which shells of 2SF phase intercalate SCF-phase shells, 
depending on the value of the local hlling factor. The different panels of Fig. 2 show density 
prohles for diherent values oi N = Na + but same parameters U/t = 20, Uab/t = 17 
and un/t = 0.12. From top to bottom, radial prohles of total density with increasing boson 
number are displayed. By increasing N, the spatial shell structure appears, revealing the 
alternate shells of the two phases, with the SCF characterized by a density plateau. We do 
not observe a plateau at n = 2 in agreement with reference [5] where the formation of SCF 
phase at double-commensurate hlling is observed for larger values of U. 
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FIG. 1. dSF phase {U/t = 10, Uab/t = 15, Na = 639, N^, = 641) confined in a Harmonic trap 
(uu/t = 0.12). Left-Panel shows density maps for both species over the entire trap. Right-Panel 
shows a radial section of density maps, perpendicular to the boundary. 
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FIG. 2. Radial profiles of total density riy for increasing number of particles at 17/t = 20, Uab/t = 17 
and oJu/t = 0.12. 
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